o 

> 
O 

o 

(N 
Oh- 



c*2 

O 
Ph. 



> 



X 



Theoretical and numerical studies of wave-packet propagation in tokamak 
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Theoretical and numerical studies of wave-packet propagation are presented to analyze the time varying 2D 
mode structures of electrostatic fluctuations in tokamak plasmas, using general flux coordinates. Instead of 
solving the 2D wave equations directly, the solution of the initial value problem is used to obtain the 2D 
mode structure, following the propagation of wave-packets generated by a source and reconstructing the time 
varying field. As application, the 2D WKB method is applied to investigate the shaping effects (elongation 
and triangularity) of tokamak geometry on the lower hybrid wave propagation and absorbtion. Meanwhile, 
the mode structure decomposition (MSD) method is used to handle the boundary conditions and simplify the 
2D problem to two nested ID problems. The MSD method is related to that discussed earlier by Zonca and 
Chen [Phys. Fluids B 5, 3668 (1993)], and reduces to the well-known "ballooning formalism" [J. W. Connor, 
R. J. Hastie, and J. B. Taylor, Phys. Rev. Lett. 40, 396 (1978)], when spatial scale separation applies. This 
method is used to investigate the time varying 2D electrostatic ITG mode structure with a mixed WKB-full- 
wave technique. The time varying field pattern is reconstructed and the time asymptotic structure of the 
wave-packet propagation gives the 2D cigenmode and the corresponding eigenvalue. As a general approach to 
investigate 2D mode structures in tokamak plasmas, our method also applies for electromagnetic waves with 
general source/sink terms, either by an internal/external antenna or nonlinear wave interaction with zonal 
structures. 



I. INTRODUCTION AND MOTIVATION 

The two dimensional mode structure (in the radial and 
poloidal directions) in tokamak plasmas and its spatio- 
temporal evolution is important since it is related to 
many physics aspects of fusion interest. For investigation 
of Alfven waves 1-6 and drift waves 7-10 , the poloidal vari- 
ation of the plasma equilibrium gives rise to "ballooning- 
like" parallel mode structures, while the instability of the 
modes, in turn, may set the maximum /? attainable in a 
tokamak 11 . On the other hand, the radial fluctuation 
structure, e.g. the radial envelope of the mode, is related 
to the radial correlation length of the underlying turbu- 
lence and might impact turbulent transport in various 
ways. Meanwhile, for radio-frequency (RF) wave heat- 
ing/current drive, the 2D mode structure, which bears 
the information of wave vector, wave intensity and loca- 
tion of the reflection/resonant /mode conversion layers, is 
crucial for calculating the heating/current drive efficiency 
and deposition profile. For the case of the lower hybrid 
wave propagation in tokamak plasmas, the toroidal ef- 
fect, i.e. the breakdown of the poloidal symmetry of the 
system, causes a shift of the parallel wave number and 
its evolution 12,13 , which is related to the current drive 
efficiency and the wave energy deposition layer. 

Besides the spatial structure, the temporal evolution 
of the fluctuations also incorporates important features. 
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For example, the Ion Temperature Gradient (ITG) mode, 
modulated by zonal flow, shows regular or chaotic be- 
haviors, depending on the relative ordering of differ- 
ent characteristic times, e.g. the inverse linear growth 
rate and the libration/rotation period of the wave-packet 
propagation 14,15 . Zonal flows are also known to be con- 
nected with ITG turbulence spreading 16,17 , which may 
occur via soliton formation 18 . These findings are help- 
ful to understand basic features of transport process, e.g. 
the transport scaling with the system size 16,17,19 . 

Although mode structures can be readily obtained 
by many numerical solvers, either by eigenvalue ap- 
proach, e.g. ERATO 20 and MARS 21 , or solving the ini- 
tial value problem as, e.g., in GTC 22 and in the flux- 
tube simulations 23 , analytical and semi-analytical tech- 
niques are developed to simplify the original problem, 
while maintaining the key physics. The ballooning for- 
malism makes use of the "translational invariance" of 
poloidal harmonics in the Fourier decomposition of the 
fluctuations, in order to reduce the 2D eigenvalue prob- 
lem to a local ID parallel eigenvalue problem; thus, it 
gives the local eigenvalue and the parallel (to B) mode 
structure 11 ' 24-28 . At the next order, the ballooning for- 
malism takes into account the mode radial envelope and 
the radial structure can be also obtained 1-10 ' 29-31 . While 
these methods arc based on the peculiar properties of 
MHD and drift wave fluctuations in strongly magnetized 
plasmas, the WKB approach reduces the problem to the 
first order partial differential equation of the eikonal, so 
that the phase and the amplitude of the wave can be 
obtained easily along the characteristic lines in weakly 



2 



non- uniform systems (see [32,33] for a review). 

In this paper, we formulate the "initial-value" problem 
of wave-packet propagation with a source term to get the 
mode structure in general tokamak geometry using flux 
coordinates. Rather than solving the partial differential 
equation for the wave directly, we can "follow" the evolu- 
tion of wave-packets generated by the source and calcu- 
late fluctuation patterns generated by them till the mode 
structure sets up. Here, the "source" can be generic; for 
example, in the case of RF wave propagation, the source 
is provided in the form of the initial phase and amplitude 
determined by the antenna located at the boundary 13,34 . 
In addition to this "external" antenna, the source term 
can represent an "internal antenna" , which can be used 
to investigate mode structure, frequency and damping 
rate as, e.g., in the case of Alfven Eigenmodes 35-37 . By 
adding the antenna induced perturbed scalar potential 
5$ ant (r,t) = A ant (r)cos(n ant 4> - m ant 8) into the origi- 
nal equation, where 4>, 8, n, m denote the toroidal and 
poloidal angles and the corresponding mode numbers, re- 
spectively, and the subscript "ant" stands for "antenna" , 
the eigenmode is excited and the saturated wave ampli- 
tude is given by 38 
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(i) 



where u) an t is the "antenna" frequency, uj 2 = w, 2 + j 2 , u r 
is the real frequency of the eigenmode and 7 is the damp- 
ing rate. With all the information provided by the an- 
tenna excitation, we can reconstruct the mode structure 
in this region. In this paper, we will discuss the gen- 
eral two dimensional propagation of wave-packets with 
an "internal antenna" inside the plasma or an "exter- 
nal antenna" at the boundary. Generally, the source can 
also account for nonlinear wave interactions with zonal 
structures 14,18 ' 39,40 and, in the framework of general co- 
ordinates, it is straightforward to extend this approach 
to three dimension when the toroidal symmetry breaks 
down. 

The paper is organized as follows. In Section II, we 
introduce the general flux coordinate system and mode 
structure decomposition (MSD) approach. There, we dis- 
cuss the general use of MSD to represent the problem in 
the mapping space and to obtain the physics solution in 
real space from that in the mapping space. The mixed 
WKB-full-wave approach for solving the 2D problem is 
also discussed along with its connection to the balloon- 
ing formalism. In Section III, we introduce the WKB 
method for studying wave-packet propagation in general 
geometry, presenting its application to the propagation of 
cold lower hybrid waves. Meanwhile, the ITG eigenmode 
formation using the mixed WKB-full-wave approach is 
analyzed and demonstrated. In Section IV, we give our 
conclusions and final discussions. Two appendices are de- 
voted to more formal discussions of the mode structure 
decomposition approach and its connection with the well 
known "ballooning formalism" . 



II. GENERAL FLUX COORDINATES AND THE MODE 
STRUCTURE DECOMPOSITION APPROACH 

In this paper, we use the general magnetic flux sur- 
face coordinates for our analyses of wave-packet prop- 
agation in plasmas of fusion interest, characterized by 
complex geometries. Some of our discussion follows R. 
White's treatment of straight field line coordinates 41 . We 
use r, 8, ( to describe the radial-like, poloidal-like and 
toroidal-like coordinates, of which the latter two have 2ir 
periodicity. The radial-like variable r is defined by 
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where a is the normalization length, e.g. the tokamak 
minor radius, ^ and arc poloidal and toroidal flux, 
respectively, while the subscripts "6" and "0" stand for 
the boundary and on-axis value, respectively In order to 
obtain the straight field line coordinates (r, 6, (), we shift 
( with respect to the toroidal coordinate <fi by v(r, 9) = 
<j) - C 41 so that 
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where v is a periodic function of 8. Here, the toroidal 
symmetry is preserved for £ since 
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according to the chain rule. The choice of 8 is used to 
choose a convenient form for the Jacobian 

J= (V r xV«.V0-' = |x|.|. (5) 



Examples are Hamada coordinates 

J = J H (r) 
and Boozer coordinates 

J=J B (r)/B 2 



(0) 



(7) 



where the Jacobian or the product of the Jacobian with 
B 2 is a function of 'J', respectively. Using equation (4), 
we can obtain straight field line coordinates where the 
magnetic field can be written as 

B=V(xV$- q(r)V8 x V* . (8) 
With the definition 

Z = (-q(r)8 , (9) 
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we can obtain the Clebsch representation for the mag- 
netic field 

B = V£xW (10) 

in the Clebsch coordinates (r, (?,£) 23 or (r, £, C) 42 i with 
£ = —£/<?, depending on the choice of the parallel co- 
ordinate along the magnetic field line as or £ and the 
rescaling of \& and £, which are stream functions of the 
magnetic field. 

Straight field line coordinates have good features, 
such as constant safety factor on magnetic flux sur- 
faces and the advantage of easily describing the paral- 
lel (to B) mode structure. In addition, different coordi- 
nate systems represent different features. For example, 
(r, £, C) can handle the geometry with q = oo (X-point 
configuration) 42 , where the (r, #,£) coordinates fail since 
£ = C — q@ ~ —q6. Usually, the appropriate coordinate 
system is chosen according to the needs of the calcula- 
tion. For example, the gyro-kinetic turbulence code GTC 
makes use of straight field line coordinates to suitably 
describe the field-aligned structure of micro-instability 22 
and the particles are pushed with guiding-center Hamil- 
tonian formalism with conserved phase volume 43 . The 
grid points in the Clebsch coordinates (r, £ , Q are aligned 
to follow the equilibrium magnetic field lines and, thus, 
the number of grid points aligned along the parallel coor- 
dinate can be greatly reduced since the modes have elon- 
gated structures along the field lines (fen <C k±). As for 
the particle pushing, a simple Rungc-Kutta scheme can 
be used safely and accurately to integrate the Hamilto- 
nian system since particles move mainly along a straight 
line in the straight field line coordinate system (r, 6,£). In 
flux-tube simulations 23 , Clebsch coordinates (r,0,£) are 
used and the global domain is replaced with flux tubes 
along the magnetic field line, * <E (*o - A*, * + A#), 
£ e (Co - A£, £ + AO, 9 e (0 - AO, 8 Q + AO). The main 
concept to deal with the boundary condition is the statis- 
tically motivated periodicity, which assumes that the sta- 
tistical properties of fluctuations at the boundaries along 
i&, £ and are identical when the simulation domain is 
larger than the correlation length along the three coor- 
dinates. When the parallel correlation length increases 
and becomes larger than 2ir, the parallel length of the 
box needs to be extended beyond 2ir and non-physical 
parallel wavelengths are permitted, which requires a ded- 
icated discussion about the global consistency 44 . Many 
works 22,23, 41-44 , including recent ones 45 , deal with the 
formulation of field-aligned flux coordinates and their im- 
plementation in numerical codes; so, we do not discuss 
this issue further. 

In our approach, in order to deal with boundary con- 
ditions, we make use of the mode structure decompo- 
sition (MSD) approach 15 , which can be reduced to the 
well-known "ballooning formalism" 11,24 ~ 30 ' 46 when spa- 
tial scale separation applies between fluctuation struc- 
tures and equilibrium non-uniformities. Writing the fluc- 



tuation in the form 

/ = e"<f n {r, 0) = e"< ]T e~ imS f n , m {r) , (11) 

m 

by continuing f n ,m{r) to tp n (r,x), i.e. <p n (r,x) € 
{(p n (r,x)\ip n (r,m) = f n ,m(r),x e K,m G Z}, and defin- 
ing the Fourier transform 

f n (r, V ) = (27T)- 1 J e-"^ l p n (r,x)dx , (12) 

we have 

f„(r,0) =2^/ rl (r,#-27rf) 

i 

= Y,e- ime I e im V„M)*7 , (13) 

m ■* 

where the Poisson summation formula has been used in 
the equivalent form of equation (Al) in Appendix A. 
Then, f n (r, 0) can be obtained by solving the linear Par- 
tial Differential Equation (PDE) 

C(r,r);d r ,d v )f n (r,r)) = , (14) 

with suitable boundary conditions in the mapping space 
and performing summation via equation (13), rather than 
solving the original equation 

C(r,0;d r ,d § )f n (r,0) = O , (15) 

with 27r periodic boundary condition in direction. In 
equations (11) to (15), the toroidal symmetry has been 
used for the reduction of the PDE corresponding to the 
considered wave equation from 3D to 2D. The partial 
differential operators above are carried out holding the 
other coordinate constant as well as the third one, e.g. £ 
in the following Clebsch coordinates. The representation 
of C(r, r)] d r , d^) can be obtained by mapping the differen- 
tial operators in C from real space to the mapping space. 
Since these operators remain the same in real and map- 
ping space (equation (A4) in Appendix A) except the 
substitution of with ?y, C remains formally the same 
in equation (14) and (15). Although the solution in real 
space can be obtained from that in the mapping space, 
in other words, equation (14) implies equation (15), the 
opposite is not true. For details, please see Appendix 
A, where equation (A13) demonstrates this statement, 
which is connected with the non-uniqueness of the con- 
struction of f n (r,0) from f n {f,9) in equation (12), for 
the sampling operator ip n (r,x) H> / n , m (r) admits dif- 
ferent choices of (p n {r,x). Some previous works 31 ' 46,47 
have discussed the uniqueness of one specific choice of 
the inverse "ballooning" transformation. However, gen- 
erally speaking, the representation of equation (12) is not 
unique and different constructions of ip n (r,x) are possi- 
ble, as discussed in Reference [15] and Appendix A. 
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In practical applications, it is often useful to adopt 
Clebsch coordinates (r, 6, £), where 

/M,£)=5> m «F„M) . (16) 

n 

Prom a direct comparison with equation (11), we obtain 
f n (r, V ) = e-"^F n (r, V ) , (17) 
from which, noting equation (13), we have 

F n (r, 9) = 2nJ2e 27rlinq F n (r, 9 - 2ir£) 

e 

= S^ e i{nq-m)6 j e ^-nq)V F n (r, V )dr] .(18) 

The original differential equation, equation (15), after 
transformation to Clebsch coordinates, can be written as 

g{r,6;d r ,d§)F n (r,6) = , (19) 

in real space, or 

g(r,md r ,d ri )F n (r,r,) = , (20) 

in the mapping space, where the linear differential oper- 
ators Q are formally the same in the two equations. The 
solution F n (r, 9) can be obtained by solving equation (20) 
with suitable boundary conditions and the construction 
by equation (18). 

The mixed WKB-full-wave approach 15,48 , as general- 
ization of the standard ballooning- mode formalism 11 , can 
be used to solve equation (20). When its applicabil- 
ity condition is satisfied, i.e. the wave propagates much 
faster in the parallel than in the radial (magnetic flux) di- 
rection, the parallel mode structure is formed before the 
wave propagates significantly away from the considered 
magnetic field line; thus, it is meaningful to separately 
calculate the parallel mode structure and the radial en- 
velope. The mixed WKB-full-wave approach, proposed 
here, is essentially the same as that used in [29-31,49,50] 
for investigating high mode number MHD modes and re- 
sulting remarkably good even down to low mode num- 
bers. The same method was used later for analyzing 
eigenmode stability of drift waves 7-10 as well as Alfvcnic 
modes 1-6 . The main novelty of our present treatment 
stands in the solution of the wave equation as initial value 
problem in the presence of a source term, which makes 
extensions to nonlinear problems readily available 15 and 
allows us to present the wave-packet propagation in mag- 
netized plasmas with general equilibrium geometries as 
one single coherent framework, with different possible ap- 
plications to propagation of radio-frequency (RF) waves 
in toroidal geometries (section III A), mode structure and 
stability of drift waves (section IIIB) as well as Alfvenic 
fluctuations and MHD modes. Since rj denotes the coor- 
dinate along the equilibrium magnetic field, as an "ex- 
tended poloidal angle" 11 , equation (20) gives the parallel 
wave equation by substituting d r with ik r (r) 

Q(r,r 1 ;ikr,d ri )A(r)F nO (r,r]) = , (21) 



where A(r) — e 1 $ krdr is written in the eikonal form. 
Then, given the solution of the parallel wave equation, 
k r is obtained from 

D(u,r,kr)A(r) = , (22) 

where ui is the mode frequency and the local dispersion 
relation D{u>,r,k r ) is obtained from the solution of equa- 
tion (21) with proper boundary conditions 49,50 . The 
radial envelope A(r) can be obtained using the WKB 
method when \d r k r /k^\ <C 1 (see section III). The appli- 
cability condition of the mixed WKB-full-wave method 
is more stringent than the obvious condition on the ra- 
tio of radial to parallel wave-packet group velocities be- 
ing small; |v g rA'g||l ~ "C 1. The condition 
that the parallel mode structure is formed before any 
significant radial propagation takes place is given by 
\v gr \ <C \v 9 \\Lt/ L\\\, with Lt ~ r the characteristic arc- 
length of the wave-packet trajectory in a toroidal plasma 
cross-section and Lit ~ qRo the connection-length of a 
tokamak of major radius Rq. Thus, the actual applica- 
bility condition of the mixed WKB-full-wave method is 
\Vgr/v g \\\ « |fc||/fcx| < r/(qR ) < 1. The mixed WKB- 
full-wave method can be used, e.g. to calculate the 2D 
structure of TAE 2 and ITG modes 7 . It has been also 
proposed to investigate the lower hybrid wave propaga- 
tion with fen <C fcj_, or more generally, with a small ratio 
of radial to parallel group velocities 48 . As an example 
of the mixed WKB-full-wave approach, used for solving 
an initial value wave equation in the presence of a source 
term, we discuss the calculation of the 2D electrostatic 
ITG mode structure in section IIIB. Meanwhile, for il- 
lustrating the calculation of the time-evolving 2D wave 
structure of RF waves in general toroidal geometry by 
means of 2D WKB method, we analyze the case of the 
propagation of a cold Lower Hybrid wave-packet in sec- 
tion III A. The application of the mixed WKB-full-wave 
approach to the LH propagation requires a dedicated dis- 
cussion and will be reported elsewhere. 

Before discussing specific applications of the mode 
structure decomposition (MSD) method, we give a brief 
discussion of its connection with the ballooning formal- 
ism. A more formal and detailed analysis is given in 
Appendix B. As anticipated above, the mode struc- 
ture decomposition method is general and can be re- 
duced to the ballooning formalism when spatial scale 
separation applies between fluctuation structures and 
equilibrium non-uniformities. In the mapping space of 
Clebsch coordinates (r, #,£), when the equilibrium pro- 
file variation is much less than that of the radial wave- 
length, then \d r F(r,r))\ -C \nq' F n {r,rf)\. This condi- 
tion is equivalent to the quasi translational invariancc 
of poloidal harmonics in ballooning formalism, where 
fn,m(r) = A(r)f n (nq + m), with A(r) being the slowly 
varying amplitude. Thus using F n (r,r]) = A(r)F n (r,r]), 
we have the ballooning representation 

f n (r,6) = A(r)J2e- im8 f e i( - m -^F n {r,r,)dr, (23) 
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where A(r) and F n (r,rj) vary on the envelope scale La 
and the profile characteristic length L p , respectively, with 
La <C L p . Using the scale separation argument, we can 
also get the inverse of ballooning transformation pro- 
posed by Hazeltine [46,47], which makes the mapping 
between f n (r,r]) and fn(f,9) to be one-to-one. In Ap- 
pendix B, we demonstrate that equation (14) is the nec- 
essary and sufficient condition for equation (15) to be 
satisfied, with proper scale separation and using Hazel- 
tine's prescription for construction of the image function 
in the mapping space. We also clarify the more general 
situation in the framework of the MSD approach. 

An advantage of moving to the mapping space to solve 
the equation for /„ (r, rf) or F n (r, rf) is that we can easily 
choose proper boundary condition in rj direction, rather 
than the 2tt periodic constraint in 6 direction. Typical 
boundary conditions are outgoing wave, where no energy 
is generated at large rj. When the "flux tube" is chosen 
long enough, the vanishing boundary condition can be 
chosen as well. Then, with the solution in the mapping 
space, the physical solution satisfying periodicity con- 
straints in 9 direction can be uniquely constructed with 
equations (13) or (18). In addition, since the MSD ap- 
proach only relies on the Poisson summation formula, it 
is not limited by the condition q' ^ required by the 
ballooning formalism, and its application to cases with 
vanishing magnetic shear is readily obtained ' . 



III. COMPLEX WKB FORMULATION OF THE 
WAVE-PACKET PROPAGATION PROBLEM IN 
GENERAL TOROIDAL GEOMETRY 

The "WKB" is a method for calculating approximate 
solutions of linear partial differential equations with vary- 
ing coefficients. In plasma wave propagation, it is applied 
to investigate the solution of the Maxwell- Vlasov system 
. Based on the eikonal form of the solution of the wave 
equation 



$(r) = Aexp{iS } , 



(24) 



it relies on the geometric optics approximation, i.e. the 
ratio between wavelength and the characteristic length of 
the variation of the reflective index must be small, 



Vk 



< 1 



(25) 



By using the ansatz equation (24) and (25), the ray 
trajectory equation system and the amplitude equation 
are obtained from the asymptotic expansion in e of the 
governing wave equation. This set of first order dif- 
ferential equations describe, at the lowest order, the 
wave propagation and, at the higher order, the focus- 
ing/defocusing effects of the propagating rays. While 
traditional geometrical optics deals with the diffrac- 
tiveless wave fields, beam tracing methods were devel- 
oped to investigate the diffraction of lower hybrid wave 



propagation 53 . The beam tracing method, even though 
it preserves the Hamiltonian particle description for the- 
wave in the group velocity v g direction as the usual 
WKB method, retains the full wave description in the 
direction perpendicular to v s ; the additional ordering be- 
tween the wavelength and wave beam size is then used 
to expand the original wave equation as asymptotic se- 
ries. The beam tracing is a complex WKB method with 
extra ordering suited to describe the finite beam width, 
which means that the physical phenomenon of diffraction 
is included in the treatment of the wave equation. The 
general complex geometric method for investigating the 
Gaussian beams in inhomogeneous media is reviewed by 
Yu. Kravtsov 54 . Here, we'll apply it to describe plasma 
waves-packet propagation and obtain the reduction to 
the traditional WKB method when the correction from 
diffraction is negligible. As anticipated in Section II, we 
will also discuss the conditions under which wave-packet 
propagation in 2D system can be studied with a mixed 
WKB-full-wave approach. 

Using WKB method, the solution of the wave equa- 
tion is constructed along the characteristic lines, parallel 
to the wave group velocity, along which the energy is 
transported. We can apply it in real space directly or in 
the mapping space, after decomposing the wave-packet 
with the mode structure decomposition (MSD) method 
(see Section II). According to the results of section II, the 
equations in both spaces are equivalent and formally the 
same. However, the different boundary conditions lead 
to different solutions connected with periodic operators 
(see Appendix A and B). 

Here, we discuss how to reconstruct the time-varying 
wave field by the WKB method. In the electrostatic limit, 
the wave equation reduces to 



£(r,V)$(r) = , 



(26) 



where $ is the perturbed scalar potential and the op- 
erator £(r, V) is determined by the plasma parameters 
and the structure of the dielectric tensor. Assuming 

V$(r) = i jv"S(r) j $(r) = ik$(r) and solving equation 

(26) with an asymptotic expansion in e, we have 

Do(r,ko)$(r)=£(r,iko)#(r) = (27) 

at the lowest order O(e ), where we have denoted with 
ko the first term of the asymptotic series k = ko + ki + 
k2 . . .. Equation (27) is readily solved by the ray tracing 
equation system (method of characteristics): 

dr _ dD a 
dr dk 
dko dD 
dr dr ' 

where t is a time-like coordinate parameterizing the 
wave-packet motion along the ray trajectory. Then, the 
lowest order eikonal in equation (24) can be obtained by 
integration along the trajectory 

r dr 

So= / ko • j^dr . (29) 



(28) 
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Meanwhile, ki, k2 . . . can be obtained by expanding the 
wave equation at the higher order. The ki correction 
corresponds to the leading order amplitude term in the 
eikonal representation of equation (24) and is given by 



A(r) = cxp{i / ki • — — g?t} 



(30) 



Knowing A and So, it is possible to reconstruct the lead- 
ing order wave field accurate to 0(e 1 ). 

In equation (28), the complex ray tracing method takes 
D as a complex function and thus, the imaginary part of 
dDg/dko leads to imaginary displacement of the ray in 
the complex space. In the framework of complex ray trac- 
ing, rays have a real physical meaning only when they in- 
tersect the real space. Thus, rays should be started with 
given boundary condition in complex space, analytically 
continued from the real boundary condition, and propa- 
gate in the complex space till they reach the real space 54 . 
Unlike beam tracing, where an additional small param- 
eter, the ratio between wavelength and beam width, is 
introduced to expand the wave equation and to force the 
beam to move in real space, the complex WKB method 
deals with the ray equation in complex space without 
additional scale separation 55,56 . Using complex WKB 
method, diffraction can be described well and a heuris- 
tic result in homogeneous media shows that the critical 
propagation length for diffraction is 



W/2 

A 



(31) 



where W is the Gaussian beam width and A is the wave- 
length. For propagation length less than L c the diffrac- 
tion effect is negligible and traditional WKB method 
works reasonably well. 

Coming to the case of tokamak plasmas, considering 
the equilibrium toroidal symmetry, linear mode struc- 
tures can be Fourier decomposed in the toroidal direction 
and, for a given toroidal mode number n, the problem is 
reduced to two dimensional and equation (26) reduces to 



£(r,9,d r ,d s )$ n (r,9) =0 



(32) 



where we assumed $(r) = e m ^$ rl (r, 6). To study this 
equation, we can generally use a 2D WKB method, as we 
will do in section III A for the case of a cold lower hybrid 
wave-packet; or we can adopt the mixed WKB-full-wave 
approach, introduced in section II, as in the application 
to the electrostatic ITG propagation discussed in section 
IIIB. In the following section III A, we will employ the 
traditional WKB method to discuss the case of the lower 
hybrid wave propagation in the electrostatic limit, recon- 
structing the time dependent wave-field pattern in 2D 
geometry. Meanwhile, to demonstrate the mixed WKB- 
full-wave method, we discuss its application in the case 
of the electrostatic ITG propagation and eigenmode for- 
mation in torus in section IIIB. The application of the 
this method to the lower hybrid wave propagation will be 



reported elsewhere, since it requires a dedicated work to 
discuss its applicability and to compare its results with 
the findings from a more conventional 2D WKB method. 



A. Propagation of a cold lower hybrid wave-packet 

As the first application of the method illustrated 
above, we consider the propagation of the cold lower hy- 
brid wave in the frequency range u> C i <C lo *C w ce . As 
a quasi-electrostatic wave, lower hybrid wave can be de- 
scribed by the Poisson equation 



£(r, V)$(r) = V • (e(r) ■ V$(r)) = 



(33) 



where e(r) = SI+ (P — S)hh — iDh ■ e is the cold dielec- 
tric tensor, e is the Levi-Civita tensor, and S,D,P are 
elements of the cold dielectric tensor in Stix notation. 
Thus we have 

V • (SV$(r)) + V • ((P - S")V||$(r)) = (34) 

and correspondingly, 

D (r,ko) = Sk* + (P - S)kf ]0 = , (35) 

Di(r,ki) = -iV-(Pk|| + Skj.o) (36) 
+2(Pk|| -k t | 1 +Sk ±0 -k ±1 ) = 

in the zeroth and first order WKB expansion. In the 
straight field line coordinates (r, mentioned in sec- 
tion II, equation (35) can be written as 



D = Sg af} k a0 kp 



P-S 



(nq + mf = , (37) 



where a,/3 <E {r, 9, £} and fc a0 , kpo are their conjugate 
momenta defined by k a0 = dSo/da, kp = dSo/d(3. The 
metric elements are defined as g a @ = e Q • with the 
contra- variant basis e Q = Vq. The repeated superscript 
and subscript mean summation. Then the ray tracing 
equations are 



^ = -2Sg r % 



(10 

dc 

(It 
dk a0 
dr 



Sg §l3 k po + (P-S) 
Sg/Vkpo + (P~S) 



nq 



J 2 B 2 
nq + m 



-q 



da ^} ' 9 da 
1 d(P - S) 



P-S 
_l_d_ 
7Bda~ 



dc 



kpokjo 
i d 



(38) 
(39) 
(40) 
(41) 



(JB) 



(P-S) 



nq -f 
(nq 



da 



J 2 B 2 



The initial condition is given at r = on the initial 
surface Q, paramctrically defined as r = Yj(l, c) where t, 
c, are curvilinear coordinates on Q and subscript / means 
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initial value. In tokamak, we can choose i = 6, <; = £ at 
the starting surface r = r(r = n,6,Q. Then the initial 
condition of cikonal Sqj and amplitude A is 



Sq\q = S i{0, C) , 
A \ Q = A 0I (9X) , 



(42) 
(43) 



and (fc r , m, ti)\q are obtained from Sqi by their definition 
and local dispersion relation. 

At the higher order, equation (36) for the amplitude 
can be written as 



dA A dr 
^ + 2 V '^=° 



(44) 



where 



dA dr 
dr dr 



VA 



dr 

dr 



dP 
dk 



To derive V • J^, we move to the ray coordinates (t, £, r), 
where t,<r label the rays starting from the initial surface 
Q and r indicates the distance along a fixed ray. In this 
coordinate, the amplitude is 



A(t)=A(t )< 



' Mtq) 
Mr) 



(45) 



where the Jacobian of the ray coordinates is defined as 



Jr. = jpp x j£ ■ §| . Applying this result to the tokamak 
geometry and making use of equation (5), we have 



Jr. = J 



dr 89 89 dr 
dr di dr dt 



(46) 



where we make use of the toroidal symmetry to simplify 
Jr since ^ = S aq . Here 5 ai is the Kronecker delta. 
This is also the formula implemented in the numerical 
calculation. 

Although the equation set is quite general, we can still 
get some qualitative information on the properties of its 
solutions, such as on the location of the reflection points 
where dr/dr = 0. From equation (38), we have 



dr 



-2Se r 



"_L0 



where e- 1 is the perpendicular unit vector. Thus two 
types of reflection points exist. One appears near the 
cut off layer with kj_o — > 0. The other appears usually 
in the inner region where e r • = 0. The lower hybrid 
wave, once it has been efficiently coupled with the plasma 
from the external launching system, is trapped between 
the two reflection points until absorbed by electron Lan- 
dau damping (ELD). This phenomenon is referred to as 
"multi-reflection" 34 . We can also assume simplified toka- 
mak geometry, e.g. concentric circular magnetic surfaces, 



to reduce the equations and obtain asymptotic solutions 
of the wave propagation semi-analytically or numerically. 
This is investigated and discussed in Ref. [57]. 

In order to calculate the propagation and absorption of 
the lower hybrid wave in a general self-consistent plasma 
equilibrium, the magnetic field B, the Jacobian J and the 
metric tensor g lJ are calculated for a given equilibrium, 
using the flux coordinate system introduced in section 
II. The determination of the region, where the wave de- 
posits its energy, is important in Lower Hybrid Current 
Drive. The central energy deposition seems to be pre- 
vented owing to the Electron Landau Damping when the 
wave propagates in high-temperature plasma or the par- 
allel wave number is large. The parallel wave number, 
fixed by the external antenna and evolving during the 
wave propagation, is a crucial parameter for the determi- 
nation of the power deposition profile. The general WKB 
formulation of the wave equation, as described above, is 
suitable for studying the effects of the plasma equilibrium 
parameters on the power deposition profile, because the 
plasma equilibrium is automatically taken into account 
in the ray equations for the phase and the field ampli- 
tude by using the flux coordinate system. Moreover, if 
an analytical magnetic equilibrium is known (e.g. the 
Solov'ev equilibrium), the ray tracing equation system 
can be numerically integrated by incorporating directly 
the analytical magnetic equilibrium, avoiding an inter- 
face that takes into account the numerical equilibrium 
(Grad-Shafranov solver) and can be a source of numer- 
ical noise for the Runge-Kutta integrator. In order to 
perform this study and to test the sensitivity of the n\\ 
evolution (and consequently the power deposition profile) 
as a function of the magnetic equilibrium, the ray tracing 
equations have been solved by changing the macroscopic 
parameters that characterize the equilibrium, i.e. the 
elongation and the triangularity, in the case of FAST 58 
plasma parameters used by A. Cardinali et al. in Ref. 
[59]. When calculating the absorption, for the sake of 
simplicity, only the linear ELD has been considered. It is 
obvious that this assumption tends to overestimate the 
absorption of LH waves in tokamak-like reactors, which 
should be analyzed within the framework of the quasi- 
linear theory. Besides adopting the analytical equilib- 
rium for studies of parametric dependences of deposition 
profiles, we have compared its results with those of a 
realistic numerical equilibrium for the ITER Scenario 2 
plasma, characterized by the same macroscopic parame- 
ters such as elongation and triangularity. This is done in 
order to investigate the numerical error when interfacing 
the ray tracing system to a numerical equilibrium and 
to demonstrate the similarity between the realistic nu- 
merical equilibrium and the Solov'ev equilibrium, which 
we have adopted as analytical model for our parametric 
studies. 

In Solov'ev equilibrium, by assuming that the pres- 
sure p(*&) and the square of the poloidal current function 
F 2 (^) arc both linearly dependent on ^, the solution of 
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the Grad-Shafranov equation is given analytically as 60 

* = 5{i[i + a-^]£ + y 2 } . (47) 

lib ^ ^0 

n2 

where ij = -jjt — 1. The subscript "6" and "0" above and 

in the following parameters such as rib and no, mean the 
boundary and on-axis value respectively. The parameters 
E and D are related to the elongation k and triangular- 
ity (5. As a result, a Solov'ev equilibrium is described 
by 6 parameters, i.e. the on-axis major radius i?o, yb, 
which is the value of R 2 /Rq — 1 at the boundary, the 
boundary poloidal flux function ^b, the on-axis poloidal 
current function Fq, E and D] summarizing, the param- 
eter list is (Ro, yb, "fb, Fq,E, D). On the other hand, the 
input macroscopic parameters are the major radius R c , 
which describes the center position of the plasma in the 
mid-plane, the minor radius a, the on-axis magnetic field 
-Bo, the on-axis safety factor go, the elongation k and the 
triangularity S, i.e. (R c ,a,Bo,qo,K,6). In the following 
calculation, we use the FAST parameters 58 , i.e. i? c =1.82 
m, a=0.64 m, q(r = a) = 4 ~ 5, density no=2xl0 14 
cm~ 3 , rib=5xl0 13 cm" 3 , temperature T e0 = T! i0 =8.5 
KcV, lower hybrid frequency /=5 Ghz and the paral- 
lel wave number at the antenna niu = 2. Bq is derived 
from the equilibrium, holding T e Q = Tio and n e o = mo 
constant for different shaped equilibrium. The values of 
K and 5 are given as parameters to investigate shaping 
effects. The fitted density profile is 

/ x _ / no , if < r < r T , 

n[r) ~ \ n [1 - a n (r - r T ) 2 ] , if r T < r < 1 , 

which becomes the usual parabolic profile in the whole r 
range when rx = 0. 

The first series of ray tracing calculation aimed at 
comparing the evolution of parallel wave number and 
amplitude considering three different elongated plasmas, 
K = 1, 1.5, 2, at fixed triangularity 5 = 0, for three differ- 
ent injection angles, a) equatorial, b) top launching {tt/2) 
and c) bottom launching (— 7r/2). Numerical results are 
illustrated in the sequence of figures (l)-(3). First, figure 
(1) demonstrates the evolution of the parallel wave num- 
ber along the ray trajectory and the absorption. The 
absorption's contribution to the amplitude is described 
by 

A abs = cxp{ / drj es } , (48) 

Jr A 

where 

7 „ = _^ exp / (^_y\ (49) 

V eth n \\ { \n\\VtheJ J 

corresponding to the given dispersion relation Do = 0, 
defined above. The results in figure (1) show that the 
elongation causes a downshift of n|| for the wave launched 
from 6a = and thus leads to more central absorption. 



However, the effects is the opposite for waves launched 
from ±7r/2, where the n\\ shifts upwards and the ab- 
sorption layer moves peripherally because of elongation. 
Here, we notice that the absorption of the wave depends 
on both T e and n\\ . If the linear absorption of the wave is 
dominated by the local value of the electron temperature, 
then the variation induced by the equilibrium quantities 
on the parallel wave number is ineffective on the localiza- 
tion of the absorption layer 61 . Here, however, since the 
T e (^) profile is relatively flat, the variation of might 
play a significant role in the localization of the absorp- 
tion layer. At the next order in the WKB expansion, 
figure (2) shows the evolution of the wave owing to the 
focusing/defocusing effects. We can deduce that the re- 
flection position and the focal point, where the amplitude 
peaks, are not the same point, as in the case of the wave 
propagation in the symmetric cylindrical geometry where 

A(r)=A(r A ) ] l^ . (50) 

During a single pass, the amplitude of the wave launched 
from 6a = tt/2 peaks before reaching the reflection point 
and the peak position moves towards the boundary when 
increasing the elongation. To investigate the geometry's 
effect on the focusing and de-focusing of the rays, a 2D 
plot of the the field amplitude in the poloidal section is 
shown in figure (3). The amplitude increases/ decreases 
where the rays converges/diverges. The geometry of the 
equilibrium in fact could change the optical properties of 
the plasma that in this case is acting like a lens. 

A second comparison regards the evolution of the par- 
allel wave-number and amplitude by considering four 
different triangularity values, S = 0, 0.1, 0.2, 0.4, at 
fixed elongation k = 1, for 9 A = 0, ±ir/2. Figure (4) 
shows that higher triangularity brings an upshift of ri|| 
for 6a = and more peripheral absorption, while the 
effects are the opposite for waves launched from 7r/2. 
The effect of the triangularity on absorption of the wave 
launched from 6a = —tt/2 is not manifest, since n| in- 
creases to a level where most of the wave is absorbed in a 
thin peripheral layer near the same position and there is 
no evidence of how the triangularity influences the evo- 
lution of no. Figure (5) and (6) give the ID and 2D plot 
of the amplitude, which shows that the wave launched 
from 6 a = tt/2 is more affected by triangularity than 
those launched from 6a = 0, — tt/2. Since the triangular- 
ity makes the propagation more peripheral for 6a = tt/2 
in this case, the focal point also moves in the same direc- 
tion. 

Finally, figure (7) shows a comparison of LHW prop- 
agation in the case of ITER Scenario 2, when us- 
ing a numerical equilibrium file and a Solov'ev equilib- 
rium characterized by the same macroscopic parameters 
(R Cl a,Bo,qo,K,S). In this case, an accurate evaluation 
of the error in integrating the ray tracing equations is 
required and compared with the analytical equilibrium 
case. To verify the numerical accuracy using the numer- 
ical equilibrium, we generate the standard EQDSK file 
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FIG. 1. The effects of elongation on n\\ and the LH wave absorption. 




for the numerical Solov'ev equilibrium using the same 
parameters (R c , a, Bo,qo, k, S) and calculate the ray tra- 
jectory, parallel wave number and the relative error. The 
results from the numerical and analytical Solov'ev equi- 
librium are in good agreement. Besides that, the prop- 
agation is calculated with the numerical ITER equilib- 



rium. In figure (7), we can see that the fitted Solov'ev 
equilibrium produces similar results as the ITER numer- 
ical equilibrium. However, some differences can be ob- 
served because the fitted Solov'ev equilibrium and the 
ITER numerical one are not rigorously equivalent, ow- 
ing to the up-down asymmetry of the ITER numerical 
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FIG. 3. The 2D plot of the LH wave amplitude in equilibria with different elongation. 



equilibrium, which is, obviously, more realistic than the 
analytical Solov'ev equilibrium. In addition, the Solov'ev 
equilibrium is characterized by p(^) and F 2 ^) that vary 
linearly with \&. This is not the case for a real ITER equi- 
librium. However, the close resemblance of the results 
for both cases shows that it is reasonable to use a fitted 
Solov'ev equilibrium for a preliminary parametric investi- 
gation of the shaping effect, to be complemented by more 
detailed investigations adopting a more realistic equilib- 
rium, after the interesting parameter range is identified. 
In addition, assuming the analytical equilibrium reduces 
the numerical error in the Rungc-Kutta integrator with 
respect to the numerical equilibrium, as well the integra- 
tion time, avoiding the interpolation of the flux function 
at each step of the calculation. In figure (7), the relative 
error along the ray trajectory is plotted in the case of the 
analytical/numerical Solov'ev equilibrium and the ITER 
numerical equilibrium. The analytical Solov'ev reduces 
the numerical error by a factor of ~ 10 2 . 

This ray-tracing analysis in a general magnetic-field 
equilibrium shows the dependence of the propagation 
and absorption of LH waves on the main parameters 



of the plasma equilibrium, namely, elongation and 
triangularity, as well as on the poloidal angle of in- 
jection (equatorial, top, bottom). The results of the 
analysis can be summarized as follows. In the case 
of a flat density profile, which is relevant for ITER 
operations, it has been shown that elongation causes 
a down-shift of nu and thus more central absorption 
of the wave energy for the equatorial injection; while, 
for top/bottom injection, nn shifts upward and the 
deposition layer is more peripheral. On the contrary, 
triangularity leads to up-shift of n\\ and more peripheral 
energy deposition for equatorial injection; while, for top 
injection, nu shifts downward and the deposition layer 
is more central. Here, we also remind the conclusion by 
A. Cardinali in Rcf. [61] that the variation induced by 
the equilibrium quantities on the parallel wave number 
might be ineffective on the localization of the absorption 
layer if the linear absorption of the wave is dominated 
by the local value of the electron temperature. For 
example, if the temperature has a large gradient near 
the boundary, then most of the energy would be de- 
posited in a narrow layer, where the temperature is high 
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enough for the electron Landau damping, regardless 
the values of elongation and triangularity. The shaping 
effect on nn and the position of energy deposition layer 
is summarized in Table I, where we can see that the 
elongation and triangularity tend to have opposite 
effects. At the next order in the WKB expansion, 



the amplitude is calculated and the 2D electric field is 
reconstructed, demonstrating the focusing/defocusing of 
the rays when propagating in a shaped equilibrium. The 
comparison of the ray trajectory and evolution of n\\ in 
analytical and numerical equilibria shows the validity of 
the investigation using the Solov'ev equilibrium model 
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FIG. 7. The ray trajectories, n\\ evolution and relative numerical error of the wave propagation in numerical ITER equilibrium 
(ITER(N)), analytical Solov'ev equilibrium (Solov'ev(A)) and numerical Solov'ev equilibrium (Solov'ev(N)). 
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Geometry 


6a 


nil shift 


Absorption layer 




-tt/2 


up-shift 


peripheral 


K 





down-shift 


central 




tt/2 


down then up shift 


mainly peripheral 




-tt/2 


up/down shift 


not evident 


5 





up-shift 


peripheral 




tt/2 


down-shift 


central 



TABLE I. The effects of elongation and triangularity on the 
evolution of and the position of the enegy deposition layer. 
8a, k, 5 stand for the launching angle, elongation and trian- 
gularity, respectively. 



and also reveals advantages and drawbacks of using the 
fitted analytical equilibrium. 



B. Electrostatic ITG propagation in tokamaks and 
eigenmode formation 

As an application of the mixed WKB-full-wave ap- 
proach, we discuss the case of electrostatic ITG propa- 
gation in tokamaks and eigenmode formation, using the 
ITG wave equation in the fluid limit for a low-pressure 
tokamak plasma, with shifted circular magnetic flux sur- 
faces, as reported in Ref. [7]. Using the MSD approach 
and spatial scale separation between equilibrium profiles 
and radial mode envelope (see Section II), we decompose 
the perturbed potential $(r, 9, t) in the Ballooning For- 
malism representation 

$M,C,t)=e mC <Kr,M) , (51) 
(t>(r,6,t) =2irA(r,t)J2$ ( i ) (0 + 2pny~ mq{S+2l " T) , 



where 5<j)(9) is the parallel mode structure in the map- 
ping space (see Section II, Appendix A and B). This 
equation is equivalent to equation (23) and 5(j){n) corre- 
sponds to F(r, n) there. Writing A{r, t) in the cikonal 
representation 

A(r,t) = e -^t+ifnq'e h dr ) ( 52 ) 
we can get the equation for the parallel structure in [7] 



Ml 



i/01 



COS ( 



- e k0 )sm9] }< 



\9 



ho) 2 ] 



= 



(53) 



Here uj u = v u /(qRo), v tl = (Ti/rrii) 1 / 2 is the ion ther- 
mal speed, i?o is the tokamak major radius, r = T e /T, 
u*ni = {T l c/eB)(k x b)-(Vnj)/nj, uj*n = {T l c/eB)(k x 
b)-(VT) /Tj, oj*pi = uj* n i + w*Ti, n-i is the thermal ion 
particle density, for which we have assumed unit electric 
charge e, Tj is their temperature in energy units, uj ci is the 
thermal ion cyclotron frequency, pi = v ti juj ci is their Lar- 
mor radius, k = — iV is the wave- vector, kg ~ (—nq)/r is 



its poloidal component, b is the unit vector aligned with 
the equilibrium magnetic field, s = (r/q)(dq/dr) is the 
magnetic shear and lud = —2k$piVu/Ro is the thermal 
ion magnetic drift frequency. The subscript in 9ko de- 
notes the lowest (zeroth) order term in the asymptotic 
series expansion for 9 k - 

Equation (53) readily demonstrates the validity of the 
mixed WKB-full-wave approach, since the perpendicu- 
lar components of the group velocity are introduced by 
the Finite Larmor Radius (FLR) effect and thermal ion 
magnetic drift. Thus, the perpendicular group velocity 
is much smaller than the parallel group velocity. As a 
result, the ITG wave propagates mainly along the field 
line and circulates in the flux surface several times with- 
out significant propagation in the perpendicular direc- 
tion. The corresponding parallel mode structure can be 
obtained by solving equation (53), with proper boundary 
conditions, i.e. outgoing wave or decaying boundary con- 
ditions. The corresponding local eigenvalue cj is derived 
from the local dispersion relation 



D (r,9 k0 ,u;)A(r,t) = 



(54) 



as a function of 9ko, w and the parameters character- 
izing the local plasma equilibrium. In other words, 
lj = uj(r,9ko) from equation (54). The WKB solution 
of the amplitude (radial envelope) can be written as 



A(r,t) = A (r)A(r,t)e~ 



(55) 



where Aq(t) = e l S n i8 k odr j g ^ ne zer0 ^ n orc fer solution 
and A(r, t) is the higher order correction, which contains 
variation on the slow spatio-temporal scales only. 

In order to obtain the higher order correction for the 
radial envelope, we reconstruct the governing differential 
equation from Do(r,9ko,u>), by substitution of ui =>■ idt 
and (9fco =>■ (—i/nq')d r , and expanding it locally along the 
characteristic of the wave-packet propagation in phase 
space as 

dD ( . d 



14.40 



duj V dt 

1 d 2 D 

2 d9 2 k0 



A 



dDp 
d9 k0 



i d 
nq' dr 



i d 
nq' dr 



.4 



m 



]A I 



nq' 



dr 



A 



where 



i d 
nq' dr 



.4 



_j_d_ 

nq' dr 

——A 

nq' dr 



(56) 
S(r,t) , 

(57) 



dt and d r on the left hand side formally act on quantities 
that follow and the source S(r, t) on the right hand side 
can represent the drive due to an "internal/external" an- 
tenna or nonlinear interactions. This is equivalent to the 
differential equation for the nonlinear system in [14,40], 
where the source term is due to the drift wave-zonal flow 
nonlinear interaction. After substituting equation (55) 
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into equation (56), and assuming that the trivial fast 
time scale variation w e 



is isolated from the source 



S = Se 



we obtain the equation for the higher order correction A 



doj dt 
1 d 



dP 1 d 
d9ko inq' dr 

1 °.A 



A 



l d 2 D 



(58) 



inq' dr \inq' dr 



A d9 



in 



inq' dr 



S(r,t)/Ao 



In the square bracket, the first term is much smaller than 
the second term in the region where the WKB approach 
applies. Far from the turning points, ignoring the first 
term in the square bracket and using the ray tracing 
equations 



dr 

rf7 = 

dOko 
~dr~ 

dt 

d7 ~ 



1 dP 
nq' 86ko 

= 1 dD o 
nq' dr 

dw ' 



(59) 



we can obtain the first order correction for A as 
A(r,t) = Ao(r)A 1 (r)I(r,t)e- i " t 



(60) 



A s \ T =p+ I drS/{iAcA x ) 
o 



where 



A x = 



A^dDo/dOkp 
^dD /d9 k0 



(61) 



and the factor dDo/dOko reflects the focusing/defocusing 
effects, which is represented by the Jacobian of the ray 
coordinates Jr in the 2D WKB method for the lower 
hybrid wave propagation (sec Section III A). Mean- 
while, A(r, t) represents the even higher order corrections 
while the source's contribution is the integral along the 
characteristics of the wave-packet propagation. When 
5 = 0, the wave-packet exhibits only the free stream- 
ing propagation, described by the phase variation, focus- 
ing/defocusing effects and higher order corrections. Near 
the reflection point where OD/dOko — > and 9ko — > 9 k oT, 
WKB breaks down and the local full wave solution can 
be obtained from the local expansion of Do(r,9ko,u) at 
r T 



.dD Q d ,dD 

^\r T ,e k0T g- t + (r-r T ) — \ 



(62) 



1 



2(nq') 2 d9\ 



G> 2 P , d^_ 

rrfikoT go 



(A A) = S(r,t) 



where rj- is the turning point position where 9 k0 = 9 k oT 
and, with S = 0, it becomes equation (23) in [48]. Then, 
WKB solution in its validity region and the local full wave 



solution at the reflection point are matched asymptoti- 
cally to give the solution in the whole radial range. After 
normalization and Laplace transform in time of (AqA), 
the local full wave equation reduces to the Airy func- 
tion equation from where one can derive the connection 
formulae [48] . In the case of an isolated mode with expo- 
nentially decaying boundary conditions outside its local 
support, we readily get the phase shift of ±7r/2 between 
the incident wave and the reflected wave at the reflec- 
tion layer. In the more general case, the wave-packet can 
still propagate and tunnel through the cutoff layer, with 
given reflection and transmission coefficients, and even- 
tually reach the boundary, henceforth bouncing back and 
forth. Application of global boundary conditions allows 
computing reflection and transmission coefficients at all 
turning points. After long enough time, the eigenmodc 
structure is eventually generated if it exists. 

The envelope tracing method discussed here, is differ- 
ent from the eigenvalue approach in Refs. [ 7,8], where 
the radial eigenfunction is obtained using the decaying 
boundary condition at r ±oo. In the envelope trac- 
ing approach, the wave-packet tunneling and reflection 
at the turning points can be taken into account and be 
described with the corresponding transmission and reflec- 
tion coefficients to be determined from global boundary 
conditions. The wave-packet, on the one hand, bounces 
forth and back in the propagation region and, on the 
other hand, can reach other regions of propagation or 
cutoff. These physics can be implemented systematically 
in a numerical solution scheme, while a direct eigenvalue 
approach requires the computation of global phase inte- 
gral quantization condition on complex phase-space con- 
tours that are dependent on the eigenvalue itself. Thus, 
the initial value approach is easier to perform numerically 
than the eigenvalue approach, except for the simple case 
of an isolated mode, where the two methods are essen- 
tially equivalent and straightforward. After the parallel 
structure and radial envelope propagation are obtained, 
the time varying 2D field <f>(r, 9) can be reconstructed 
from the solution of equation (53), 5(f)(6), and the ampli- 
tude A(r, t), using equation (51). 

The 2D structure of ITG mode is important for un- 
derstanding anomalous transport. The radial structure 
of the ITG mode has been investigated in previous lit- 
erature [7-10]. Here, we go further and analyze the 
time varying 2D ITG mode structure by means of an 
initial value approach, which also makes it possible to 
compare our findings with other 2D ITG mode solvers. 
The shaping effects from the equilibrium geometry can 
also be readily investigated adopting general coordinates, 
as shown in the case of lower hybrid wave propagation 
discussed in Section III A. 

In order to calculate the 8(f) numerically, we reduce 
equation (53) to the dimensionless form 



n D n [cos 9 + s(9-9 k0 ) sin 9} 8<j> = 



(63) 
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where ft = ui/ujti, = u D /uj ti = -2kepiq, f2*„j = 
^>*ni/^u — ~kePiqRo/ L n and Q* P i = io mi /io ti = 
—kgpiqRo/Lp. We consider typical parameters at the ref- 
erence radial position ro to be s = 1, q = 2, kgpi = 0.3, 
t = 1, Rq/Lt = 3, Ro/L n = 0. The ion temper- 
ature gradient profile is assumed in the form g(x) = 
4(i?o/ Lx)(e x + e~ x )~ 2 , where x is the normalized shifted 
radial coordinates x = (r — ro)/Ar. In the following, we 
choose ro = 0.5a, Ar = 0.2a, while the toroidal mode 
number is assumed to be n = 38 from kgpi = nqpi/r rj 
2nqpi/a, where kgpi s» 0.3 and a/p^ rj 500. The safety 
factor gradient g' can be estimated as q' = sq/r « 4/o 
and thus the eikonal J drnq'6 k n ~ 30.4 J dx9ko- 

Figure (8) gives the parallel structure in the mapping 
space at x = 0. The parallel mode structure is localized 
and peaked near 6 = 6 k n as expected. Figure (9) shows 
the ion temperature gradient profile (left) and the local 
eigenvalue variation with respect to Rq/Lt (right). The 
local growth rate, in the form of the imaginary part of 
n(r,dico), increases when the ion temperature gradient 
Rq/Lt increases. The local growth rate is also affected 
by Oka, since the parallel mode structure is peaked near 
9 — 6ko ~ and, thus, for 8 k n = the mode is localized 
in the bad curvature region and is most unstable. 

When solving the radial envelope structure as initial 
value problem, the constant eigenfrequency curves gen- 
erate two types of trajectories in phase space (9 k n,r), 

I 



classified as "lib-rations" or "rotations" on the basis of 
topology arguments and investigated analytically and nu- 
merically in [2,7-10,50]. The lowest order solution for the 
normalized amplitude An (r) in the region between turn- 
ing points is shown in figure (10). The wave launched at 
the left turning point tti propagates to the right turn- 
ing point tti (left frame in figure (10)), where the phase 
is shifted by 7r/2, and then propagates leftward (center 
frame in figure (10)). The interference patterns are gen- 
erated because the wave bounces between the turning 
points, which, in this case, are located at x = ±0.3 for 
a phase space "libration" (right frame in figure (10)). In 
figure (11), the next order correction for the amplitude 
A\ is plotted. As expected, its radial structure is fiat in 
most part of the propagation region around x = 0, while 
near the turning points, where dDo/d6ko — 0, the wave 
group velocity vanishes, leading to the increase of the am- 
plitude. For the sake of simplicity, only the wave-packet 
propagation region is plotted, while the exponential de- 
cay outside the two turning points and the smooth con- 
nection of the solution of equation (62) near the turning 
points is implicitly assumed. 

To illustrate the 2D time dependent ITG mode struc- 
ture, we assume a point-like source, S(r,t) = S(r — ta), 
where is the antenna position at Vti, r& = Tti- Then 
equation (60), with the phase shift at the reflection layer, 
can be rewritten as 



e wt A{r, t) =exp jz jf nq' (0 M + 6 kl ) dr - i ^ 5 ;7 r/2 j A(r, t) x (64) 

i r {t) r ^ 

As\r=o + - 5(r~r A )exp -i nq'(0 ko + 9 kl )dr + i} y 8nr/2 

1 •'0 J° i=i 



dr 



where Nt is the number of turning points the wave- 
packet passes by in the time interval t £ (0, t) and 
Si = ±1 depending on the phase shift of the reflected 
wave-packet at the turning point. Furthermore, 



6 



fci 



1 ln d Dp 

2 n d9 k0 



When propagating between the turning point pair, the 
wave-packet passes ta at t = tAo, tAi ■ ■ ■ tAN T (see fig- 
ure (12)). Assuming As\ T =o = 0, Nt as a positive odd 
integer and defining the rotation/libration number (de- 
creased by one) as N — (Nt — l)/2, we obtain 



e luJt A(r,t) =exp| l J nq'(d k0 + 8 k 



1=1 



Nt — 



i I nq (& k0 + 




i)dr-iJ2^/^jMr,t) : 

9 kl )dr + iU{N T - l)^ <W2 

kl )dr - M(N T - N T - 1) ^ <5;tt/2 



Nt—0 



l=N T + l 



(65) 
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FIG. 9. The ion temperature gradient profile (left) and the ITG local eigenvalue f2 as a function of Rq/Lt (right). 



where the Heaviside function H is used to eliminate the 
summation if the lower bound of its running index is 
larger than its upper bound. This equation reduces to 



e iut A(r,t) = A(r,t)S N (t) 



where 



i M)= ^Miexp 



nq'(8 k0 + 9 kl )dr 



(66) 



(67) 



+ cxp 



i nq'(6 k0 + 6 k i)dr - S Nt it/2 

J tA(N T -l) 



Sn= — — , e = exp{i$ } 



1 - e 

$0 = i / nq'Okodr 
= i f nq'9 kQ dr - /?7r , 



(68) 



where (3 is the phase shift (Maslov index), with /3 = 1 
for librations and /3 = for rotations. Equation (66) 



agrees with the analysis in [15] and illustrates the eigen- 
mode formation and phase mixing. In fact, Si\r(t), as a 
function of D, and N(t), describes interference patterns 
that becomes narrower around $o = 2/7T for increasing 
N, corresponding the global (2D) eigenfrequencies fie?; > 
with I € Z as the radial mode number (the left frame in 
figure (13)). The right frame in figure (13) shows that, 
unless n = fid > i-e. $0 = 2?7r, the normalized amplitude 
will decay because of phase mixing when N increases. 
The 2D mode structure with fl = (0.8582, 1.235) is se- 
lected as a eigenmode according to Sjv's time (N) evolu- 
tion. Figure (14) illustrates the 2D time dependent ITG 
mode structure. The left top frame shows the 2D mode 
structure with Q = (0.8582, 1.235) and N = (for one 
period of the oscillation in phase space) . The wave pat- 
tern and intensity remain the same for N = 8 (the top 
center frame in figure 14), since f2 = (0.8582, 1.235) is the 
cigcnfrcquency of the 2D problem. The phase mixing of 
the 2D structure with Q = (0.8576,1.231) is illustrated 
in the bottom figures. The field intensity decreases to a 
low level due to phase mixing, while the mode structure 
remains the same from N = (left bottom) to N = 8 
(center bottom) as a consequence of our present assump- 
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FIG. 10. Ao(r)'s evolution versus x. Left: propagation from the left to the right; center: propagation from the right to the 
left; left: the superposition of the previous two. 
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FIG. 11. ^4i(r)'s evolution versus x. 



tion of a point-like source. For a more general source 
term, e.g. a broad internal source, the mode structure 
will also change and short wave number structures will 
be generated due to phase mixing: the radial eigenmode 
structure is preserved asymptotically in time only for the 
correct global eigenvalues. The third column of figure 
(14) shows the fine structure for 8 £ [— 7r/6, 7r/6], where 
we can recognize the radial envelope and the variation of 
poloidal harmonics. 

In the example above, the mixed WKB-full-wave solu- 
tion for the 2D ITG mode structure is derived within the 
framework of the mode structure decomposition (MSD) 
method, which coincides with the Ballooning Formalism 
in this specific case, since spatial scale separation ap- 
plies. The parallel mode structure is obtained by solving 
the local eigenvalue problem on a given flux surface. The 
radial propagation of the envelope wave-packet has been 
solved, including focusing/defocusing effects, using the 
WKB formulae for wave-packet propagation. At last, 
the time dependent 2D mode structure is reconstructed, 
which illustrates the phase mixing and the 2D eigenmode 
formation, and also provides a picture for comparison 



FIG. 12. The radial propagation of the ITG wave-packet be- 
tween the turning points. The wave-packet is launched at ta 
and propagates between m and tti- 



with other 2D numerical solvers. A more realistic inves- 
tigation, with geometry effects for the 2D time depen- 
dent ITG mode structure and multiple space-time scale 
analysis, including the source/sink term and nonlinear in- 
teractions, can be straightforwardly implemented in the 
framework of the mixed WKB-full-wave formalism de- 
scribed here and assuming the general coordinate repre- 
sentation and MSD approach, described in Section II and 
Appendixes A and B. 



IV. CONCLUSIONS AND DISCUSSIONS 

The aim of this work was providing a theoretical frame- 
work for investigating electrostatic wave-packet propa- 
gation in tokamak plasmas with general geometry. Dif- 
ferent techniques have been discussed, ranging from the 
2D WKB approach to the mode structure decomposition 
(MSD) method and the mixed WKB-full-wave analysis, 
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FIG. 13. Sat's evolution versus Qr (left) and N (right). 




(R-R )/a (R-R )/a (R-R )/a 

FIG. 14. The time evolution of the 2D field, Real{e™V(r, 6 , t)}. Left top: ft = (0.8582,1.235), N = 0; center top: ft = 
(0.8582, 1.235), N = 8; left bottom: ft = (0.8576, 1.231), N = 0; center bottom: ft = (0.8576, 1.231), N = 8. The right figure 
show the fine mode structures for £ [n/6,n/6] in both cases. 
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which can be used for investigating RF wave propagation 
and absorption as well as the the time evolving structures 
of MHD fluctuations and drift waves. 

As application, we adopted the 2D WKB description 
to investigate equilibrium shaping effects, i.e. elongation 
and triangularity, on the lower hybrid wave propagation 
and absorbtion in tokamaks, using the Solov'ev equilib- 
rium. Numerical results show that, for midplane wave 
launching, increasing elongation leads to down shift of 
the parallel wave number and central absorbtion, while 
for top and bottom launching, the effect tends to be the 
opposite. Triangularity has the opposite effect of elonga- 
tion, except for bottom launching, where the consequence 
of increasing triangularity is not evident. The 2D mode 
structure of the lower hybrid wave is reconstructed by 
following the wave-packet propagation, including focus- 
ing/defocusing effects. These studies have been extended 
to general geometries and equilibria given numerically, 
using a reference ITER equilibrium and the correspond- 
ing "fitted" Solov'ev model. Good agreement was ob- 
tained for these cases, showing that it is reasonable to 
use a Solov'ev parameterization for identifying interest- 
ing equilibrium parameter ranges to be investigated in 
more detail with full numerical equilibria. 

As further application, we used the mixed WKB-full- 
wave method to explore the time-dependent 2D electro- 
static ITG mode structure in the presence of a (point- 
like) source term. While the parallel wave equation is cal- 
culated as an eigenvalue problem, the radial propagation 
is investigated using WKB for solving the initial value 
radial envelope problem. The 2D global eigenmode is ob- 
tained by observing the time-asymptotic behavior of the 
wave-packet, propagating between WKB turning points, 
while generic mode structures decay as the time increases 
because of phase mixing. For the sake of simplicity, we 
assumed electrostatic fluctuations for both lower hybrid 
wave and ITG mode. However, the present wave-packet 
tracing method is suitable for analyses of electromagnetic 
waves as well, including nonlinear wave interaction with 
zonal structures in the presence of a general source and 
sink. This perspective opens up the possibility of future 
applications of this present theoretical framework to a 
variety of interesting topics. 
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Appendix A: Derivation of the mode structure 
decomposition approach. 

In this appendix and in appendix B, we discuss the 
detailed aspects of the mode structure decomposition ap- 
proach, introduced in section II, while its connection with 
the well-known "ballooning formalism" (BF) 11 ' 24-30 ' 46 is 
analyzed in appendix B. 

As pointed out in section II, the mode structure de- 
composition (MSD) approach was introduced in Refs. [ 
15,48] for general analyses of wave/packet propagation in 
toroidal systems. Its validity, however, is more general 
and can be extended to any three-dimensional systems 
with two periodic dependences, one of which is ignor- 
able, i.e. corresponds to a symmetry of the system itself. 
Fusion plasmas magnetically confined in helical systems 
and planetary magnetospheres are two obvious examples. 
With respect to Refs. [15,48], we provide here a more de- 
tailed and rigorous formulation of the MSD approach and 
the derivation of its properties. Following notations in- 
troduced in section II, we use a straight magnetic field 
line toroidal coordinate system (r, 9, ^23,31,44,62,63^ wnere 
r is a radial-like coordinate depending only on the mag- 
netic flux function ip, 9 and £ are periodic angle like 
coordinates, chosen such that B • VC/B • V9 = q(r) 
and Q corresponding to the symmetry of the system, the 
equilibrium magnetic field has the Clebsch representation 
B = V£ x \7ip and £ = £ — q{r)9 64 . 

The formulation of the MSD approach and the deriva- 
tion of its properties are based on the Poisson summation 
formula, which, in its most common form, reads as 

J2e lm6 = 2tt^(S((9 -2ttto) . (Al) 

m m 

Here and in the following, the summation is implicitly 
assumed to be on m G Z. In general, using the symmetry 
properties of the system, the £ dependence of a scalar 
function /(r, 9, C), describing a generic fluctuating field, 
is decomposed as Fourier series 

/M,c) = 5> mC /nM) , (A2) 

where the Fourier components f n (r,9) are independent 
in the linear limit. Note that time dependences are as- 
sumed implicitly for simplified notations. The formula- 
tion of the MSD approach is based on the properties of 
sampling and pcriodization operators, which are related 
by the Fourier transform due to the Poisson summation 
formula, which holds pointwise in the Schwartz space, but 
holds as well in a more general context under weaker con- 
ditions 65 . The pcriodization operator /„(r, 77) \-> f n (r,9) 
from L : (R) onto L 4 (T) is expressed as 

f n {r,6) =2nY,fn(r,9-27r£) 

i 

= Y^e- me [ e im «f n (r,v)dr) , (A3) 

m J 
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where the second line is readily derived by use of equation 
(Al) and the rj integration is on R. The periodization 
operator sum converges absolutely a.e. on R 65 , which is 
a useful property for the manipulations made hereafter. 
Note that the possibility of introducing the periodization 
operator of equation (A3) was adopted in Refs. [28] for 
discussing the properties of the BF. Using equation (A3), 
it is readily recognized that 



de\ rX fn(r,9) i-> d n \ r( f n {r,r]) 



(A4) 



where p(9) is a generic periodic function. Once / n (r, vj) is 
known, f n (r, 9) is uniquely determined by equation (A3). 
The opposite is not true. In fact, the Fourier transform 
of the periodization operator f n {r,rj) H> f n (r,9) is the 
sampling operator ip n (r,x) i— ¥ (ip n (r,m)) 



with 



¥>n(r,x) = / e lvx f n {r,r))dr] 



(A5) 



Equation (A3) is readily rewritten as 
/n(M) = ^e- <m Vn(r,m) 

thus, the unique construction of 

/„(r,7 ] ) = (2 7 r)- 1 / e-^^ n (r,x)dx 



(A6) 



(A7) 



is prevented from the fact we know the function <p„(r, x) 
via the sampling operator tp n (r, x) i-> (¥>n(»", "i)) m6 z- 
This fact was noted in Refs. [15,48] and is not a problem, 
for the physical field f n (r, 9) can be always uniquely con- 
structed from / n (r, 77), while it is not necessary and often 
not even useful to construct a unique form of fn{r,rj), 
when the physical solution is known already. This view- 
point is identical to that adopted by Dewar and cowork- 
ers on the construction of the inverse ballooning trans- 
formation 31 (see appendix B). Possible constructions of 
ip n {r,x) and, therefore, of f n (r, rj) by equation (A7), are 
obtained introducing the window function w(x) 47 , which 
can be defined as a piecewise continuous function with 
maximum w(0) = 1 at x = and vanishing everywhere 
outside the interval (—1,1). In fact, we have 



V?„(r, x) = ^2 w ( x - m)ip n (r, m) 



(A8) 



/„(r, rj) = J2 e- imv 9n(r, m)— \ er^~^w(x - m)dx = g( V )f n (r, r,) 



(A9) 



where g(rj) indicates the Fourier transform of the window 
function; i.e. 

g(rj\ = — I e- lvx w(x)dx . (A10) 

27T J 

Explicit examples of the functions w(x) and g{rj) are 
given in Refs. [15,48], with g(rj) satisfying the condition 

2^]T.g(?7-2^m) = l , (All) 

which is straightforward consequence of its definition. 
I 



The periodization operator of equation (A3) and 
its non-unique inverse of equation (A9) arc of par- 
ticular interest when connected with the action 
of a generic bounded linear differential operator 
£(r,9;d r ,de)f n (r,9), which, considering equations (A4) 
and invoking the system periodicity in 9, readily maps to 
C(r,Tf',d r ,dfj)f n (r,rj). In fact, using the absolute conver- 
gence a.e. in R of the periodization operator in equation 
(A3), we have by definition: 



C(r,9;d r ,d e )f n (r,9) = / C(r,ri;d r ,dr,)f n (r,r])S(ri - 6 + 2irl)dr] . (A12) 



This equation shows that, if f n {r,rj) satisfies C(r, 7/; d r , 9 r) )/ n (r, 77) = 0, then f n {r,9) satisfies 
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£(r,9;d r ,dg)f n (r,9) = 0. Note that it was pointed out 
already in Refs. [28,31], discussing the properties of the 
BF, that f n (r,9) obeys the same equation as f n {r,rf). 
Thus, knowing the solution of a PDE in the (r, rj) space, 
allows us to construct uniquely the solution of the 
corresponding problem in the (r, 9) space. Section II 
provides a discussion of the practical advantages coming 



from this property. 

The opposite of this argument does not hold. i.e. 
given the solution of a PDE in the (r, 9) space, we can- 
not construct the solution of the corresponding prob- 
lem in the (r, 77) space, even admitting that the mapping 
/ n (V, 9) 1 — y fn(r,r]) of equation (A9) is not unique. With 
f n (r,i]) = g(ri)fn(r,i]), equation (A12) becomes 



I 

£(r,6;d r ,do)f n (r,6) = 2tt^ f £(r,ri;d r ,d n )f n (r,r])S(v - 6 + 2izl)dr, 

= 27rJ2 C ( r > r r,d r ,d e )[fn(r,6)g(9-2irt)] . (A13) 



Therefore, that £(r,9;d r ,dg)f n (r,9) = docs not 
ensure that £(r,r);d r ,d v )f n (r,r]) = but that 
X^£(r, 77; d r ,dg) [/„(V, 9)g{9 — 2tt£)] = 0, due to the con- 
dition of equation (All). Requesting that f n {r,r]) obeys 
the same equation as /„ (r, 9) yields the paradox that 
g{rj) = 1 , i.e. w(x) — 2ir6(x), thus violating equa- 
tion (All). As mentioned above and noted in previous 
literature on the BF 31 , this is not a problem, for the only 
relevant fact is that the physical field /„ (r, 9) can be al- 
ways uniquely constructed from /„ (r, 77) as solution of the 
corresponding PDE (see also appendix B). 

In practical applications, it is often useful to move from 
straight magnetic field line toroidal coordinates (r, 9, () 
to Clebsch coordinates (r, 6>, ^23,31,44,62,63^ with £ = £ _ 

q(r)9 6i . It is readily shown that equation (A2) becomes 

/M,£) = 5> m?F «M) • ( A14 ) 

n 

so that, while periodicity in £ is maintained, periodicity 
in 9 is substituted by F n (r,9 + 2tt) = e 2mnq F n (r,9). In 
(r, 77) space, the transform corresponding to using Cleb- 
sch coordinates is obtained by letting 

Ur,V)=e- inm F n {r,r,) , (A15) 



corresponding to the periodization operator 

F n (r, 9) = 2tt J2 e 2Mnq F n (r, 9 - 2itl) 

i 

= J2 e i{nq-m)e f e i(m-nq) V ^ ^ ^ 
m J 

Considering the chain rules 

d r \ g ^ = d r \ ex + q'(r)9 d c \ r 9 

de\r t £ = de\ rX +q(r) d c \ r f) 

di\ r ,e = d <\r,e > ( Al7 ) 

and that, due to the definition in equation (A15), the 
following mappings apply in the (r, 77) space 

drfn{r,r]) H> (d r - inq' (r)rj)F n (r, 77) 

d v fn(r,v) H> {d v - inq{r))F n {r,rj) , (A18) 
it is evident that £(r,9;d r ,dg)f n (r,9) M> 
Q(r,9;d ri dg)F n {r,9) and that F n {r,9) obeys the 
same equation as F n (r, 77) 28 ' 31 (cf. the above discussion 
for f n (r,0) and f n {r,rj)). In fact 



£(r, 9- d r ,dg)f n (r, 9) = e~ inq6 '£(r, 9; 8 r - mq'9, dg - inq)F n {r, 9) 
= e- mqe G(r,9;d ri dg)F n (r,9) 

= 27re~ mqd I e 2 * Unq £{r, 77; d r - inq'r], d v - inq)F n (r, r,)8{r] -9 + 2n£)dr) 
1 J 

= 2ne- inq9 Y, [ e 2Mnq g(r, m d r ,d n )F n {r, 77)^(77 - 9 + 2^dq . (A19) 
1 ^ 



As concluded above for the functions f n (r, 9) and /„(r, 77), solving the PDE 



g{r,-q;d r ,drj)F n {r,7i) = 
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and constructing F n (r,9) from equation (A16) ensures 
that Q(r, 9; d r ,dg)F n (r, 9) = 0. The opposite is not true; 
however, this fact does not pose issues of physical consis- 
tency 15 ^ 31 ' 48 . 



Appendix B: From mode structure decomposition to 
ballooning transform 

Equation (A16), which reflects the periodic structures 
underlying a generic fluctuation in a double periodic sys- 
tem that is invariant under rotations in one of the two 
periodic coordinates, is the form used in the original 
works on the ballooning representation 11 ' 24-30,46 for an- 
alyzing its properties and elucidating its usefulness in 
stability analyses of short wavelength (high-n) modes. 
In fact, when spatial scale separation between equilib- 
rium profiles and radial wavelength applies, such that 
\d r F n (r, rj)\ <?C \nq'F n (r,r])\, then the radial structures of 
e~ tnq F n (r,9) Fourier components essentially depend on 
(nq — m), thus they are characterized by a near trans- 
lational invariance. Actually, it has been noted that the 
ballooning formalism (BF) is a convenient formulation for 
treating the "flute-like" structures (\nq — m\ <C 1) that 
naturally appear as fine radial scales (inertial/resistive 
layers) in resistive 66 and ideal 67 MHD treatments of ar- 
bitrary mode number fluctuations. For the same reason, 
the same approach has been adopted for analyzing in gen- 
eral the multiple spatial-scale nature of kinetic MHD and 
Alfvenic modes excited by supra-thermal particle popu- 
lations 68-71 . The connection of equation (A16) with the 
BF and its possible generalizations has also been explored 
in situations where the local magnetic shear, rq'(r)/q(r), 
is vanishing, i.e. where the usual formulation of the BF 
does not apply 51,52 . 

The strength of the mode structure decomposition 
(MSD) approach, discussed in section II and appendix A, 
is that it introduces in general the same formal proper- 
ties of the BF, without any request of spatial scale sep- 
aration or finite magnetic shear, and it reduces readily 
to the usual BF in the appropriate limits where spatial 
scale separation applies 15,48 . Applications of the MSD 
approach are given in sections III A and IIIB. 

Equation (A16), with spatial scale separation, can 
be considered as definition of the ballooning transform. 
From previous literature and appendix A, we know that 



solving the PDE 

g(r 7 n;d r7 d ri )F n (r,n) = 

and constructing F n (r,9) from equation (A16) ensures 
that G(r,9;d r ,dg)F n (r,9) = 0. Appendix A also dis- 
cussed why, given F n (r, 9), is generally not possible to 
uniquely construct F n (r, rj) . Although this fact does not 
pose issues of physical consistency, the question of the 
uniqueness of the inverse ballooning transform has at- 
tracted significant interest 31,46,47 . The crucial ingredi- 
ent for the construction of the inverse ballooning trans- 
form 46,47 is the (infinite) separation of scales, which 
"identifies" radial coordinate and parallel wave-vector, 
i.e. introduces the notion that radial and parallel coordi- 
nates are Fourier transform duals or conjugate variables. 
In general, when spatial scale separation does not apply, 
the function ip n (r,x), entering the inverse transform def- 
inition of equation (A7), is known only via the sampling 
operator tp n (r,x) ^ {ip n (r,m)) m&z . 

Following Hazeltine and Newcomb 47 , the ballooning 
transform and its inverse are given by [cf. equation(A16)] 

f n (r,9) = f n (x,9;r) (Bl) 

= 2ir^e mq ^-^F n {r,9-2nt) , 

t 

F n (r,9) = ±- J g{s)e i ^ e f n {x + s,9;r)ds . (B2) 

Note that, here, spatial scale separation has been indi- 
cated explicitly by the notation f n (r,9) = f n (x,9;r), 
with x = nq. Furthermore, the function g is the Fourier 
transform, given in equation (A10), of the window func- 
tion w defined in appendix A. In the work by Hazeltine 
and Newcomb 47 , the function g(s) = (7rs) _1 sin(7rs). In 
the following, we show that the inverse transform, given 
in equation (B2), applies for any generic form of g(s) 
given in equation (A10). 

Since (infinite) spatial scale separation applies, we have 
that e M f n (nq + £,9;r) = f n (nq,9;r). In fact, 

f n (nq +£,0;r) = Y^ e~ ime f nm {nq + £-m;r) 

m 

= e- U8 J2e- ike f nm (nq-k;r) (B3) 

k 

= e- ue Mnq,9;r) . 

As in Ref. [47], we substitute equation (B2) into equa- 
tion (B2) and show that we have an identity: 
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f n (x,6;r) =J2e mq(2ne - 6) [ g(s)e^ + ^- 2 ^ f n (x + s,9;r)ds 
= I 9 {s)e lse Y,5{s - m)f n {x + S ,9-r)ds 
= / 9( s )^2$( s - m)f n (x,0;r)ds 
= few(2irt)j f n (x,6;r) = f n (x,6;r) . 



(B4) 



Here, we have used equation (B4) in the second line and, 
in the fourth line, the fact that J^e w(2ir£) — w(0) = 1 
by definition of the window function w, given in ap- 
pendix A. Note that, with the special choice of g(s) = 
(tts)^ 1 sin(7rs), made in Ref. [47], the identity is proved 
in the second line, since J2 m d( s )H s ~ m ) = $( s )- An 
identity is also found when substituting equation (B2) 
into equation (B2). 

K(r,0) =EJ g{sy 2 ^ x+ ^F n {r 1 8-2^)ds 

= E,e w (- 27ri ) ei2ntx ^n{r,0 - 2tt£) = F n (r,6) .(B5) 

Thus, the results reported here generalize those given 
by Hazeltine and Newcomb 47 and show that the most 
general pair of ballooning transformation and its inverse 
are provided by equations (B2) and (B2). 

The general discussion of the connection of the 
MSD approach with the standard BF is completed 
by the proof that, if f n (r,6) = e~ lnqe F n (r,8) sat- 
isfies £(r,8;d r ,de)fn(r,9) = 0, then f n (r, r)) satisfies 
£(r,r);d r ,d v )f n (r,r)) = 0. The opposite is always true 
even in the more general MSD approach, as shown by 
equation (A12). Given the inverse transform of equa- 
tion (B2), the following mappings are readily shown: 



{do 



d r f„(r, 6) ^ d r f n (x + s, 0; r) 
ix)f n (r,6) ^ (i(s + x)+d 8 )f n (x + s,9;r) (B6) 



where the second line represents the mapping of the par- 
allel derivative. The mapping of the poloidal derivative 
is trivially obtained from the spatial scale separation ar- 
gument, for which dgf n (x, 9; r) ~ —ixf n (x, 9; r). Thus 

def n {r, 9) ^ (is + d e ) f n (x + s, 9- r) (B7) 
~ —ixf n (x + s, 6; r) . 

Equations (B6) and (B8) are the formal proof that, 
if f n (nq,6;r) and F n {r,6) = e m i 9 f n (r,9) arc re- 
lated by the transforms of equations (B2) and (B2), 
C(r,6;d r ,dg)f n (r,6) = implies C{r,-q;d r ,d n )f n (r,-q) 
0. 
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